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Abstract 

We study the propagation properties of a single vortex in square Josephson- 
junction arrays (JJA) with free boundaries and subjected to an applied dc 
current. We model the dynamics of the JJA by the resistively and capacitively 
shunted junction (RCSJ) equations. For zero Stewart-McCumber parameter 
j3 c we find that the vortex always escapes from the array when it gets to the 
boundary. For (3 C > 2.5 and for low currents we find that the vortex escapes, 
while for larger currents the vortex is reflected as an antivortex at one edge 
and the antivortex as a vortex at the other, leading to a stationary vortex 
oscillatory state and to a non-zero time-averaged voltage. The escape and the 
reflection of a vortex at the array edges is qualitatively explained in terms of 
a coarse-grained model of a vortex interacting logarithmically with its image. 
For (3 C > 50 we find that the reflection regime is split up in two disconnected 
regimes separated by a second vortex escape regime. When considering an 
explicit vortex- antivortex pair in an array with periodic boundaries, we find a 
soliton-like non-destructive collision in virtually the same current regimes as 
where we find reflection of a single vortex at a free boundary; outside these 
current regimes the pair annihilates. We also discuss the case when the free 
boundaries are at 45 degrees with respect to the current direction, and thus 
the angle of incidence of the vortex to the boundaries is 45 degrees. Finally, we 
study the effect of self-induced magnetic fields (for penetration depths ranging 
from 10 to 0.3 times the lattice spacing) by taking into account the full-range 
inductance matrix of the array and find qualitatively equivalent results. We 
also discuss possible consequences of these results to experimental systems. 

PACS numbers: 74.50.+r, 74.60.Ge, 74.60.Jg, 74.70.Mq 
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I. INTRODUCTION 



In recent years the viscous motion of a single vortex in two-dimensional Josephson- 
junction arrays ( J JA) has been studied numerically by several authors |l]-[7| . In these studies 
the I-V characteristics were calculated for JJA with periodic boundary conditions perpen- 
dicular to the direction of the applied dc current. Apart from the vortex viscosity, the nature 
of its predicted mass [§-0] has attracted significant interest, also experimentally [15,16 



Van der Zant et al. reported experimental evidence for ballistic motion of vortices in a 



current-free region in a highly underdamped Josephson-j unction array [15 



An interesting aspect of this problem, which is relevant to experiments such as those in 
T5| , but that has not yet been studied systematically, is the influence of boundaries on the 



vortex motion. Experimentally one may distinguish two types of boundaries (see e.g. Ref. 
fl5fl): either the junctions at the edge are connected to a superconducting busbar or they are 
not (free boundary). The influence of the boundary on a vortex due to the busbar geometry 
can be described by an image vortex of the same sign, equidistant on the other side of the 
boundary; in the free boundary case the image of the vortex is of opposite sign. The latter 
situation is to some extent analogous to that of a soliton in a continuous long Josephson 
junction (LJJ) [17]-[19 . 



The goal of this paper is to study the reflection, transmission and annihilation properties 
of vortices in JJA with free boundaries. Here, we report simulations of dc biased arrays with 
free boundaries, in zero applied magnetic field. For every bias current considered we use 
the same initial phase configuration with one vortex in the middle of the array. We study 
the dynamical interaction of the vortex with the free boundary as it moves towards it for 
bias currents above the depinning threshold. The vortex-boundary interaction is equivalent 
to the effect of an image antivortex at equal distance on the other side of the boundary. 
A vortex and an antivortex interact logarithmically (at sufficiently large distances), so the 
same holds for a vortex and a free boundary. The total effective vortex potential is the sum 
of the interaction potential and the periodic lattice pinning potential. The Lorentz force on 
the vortex due to the applied bias current corresponds to a tilt of this effective potential. 
Below we will see this more explicitly when we discuss the coarse-grained model equations 
for the vortex dynamics. 

The outline of this paper is the following. In section II we define the model equations 
studied in this paper. There we consider the extreme type-II limit (infinite penetration 
depth) as well as the case of finite penetration depth. We also discuss the coarse-grained 
vortex equations which are used to analyze the results later in the paper. In section III (a) 
we present results in the extreme type-II regime and in III (b) the results for finite magnetic 
penetration depth. In section IV we analyze the results from III (a) in terms of the coarse- 
grained, phenomenological vortex equations including the interaction of the vortex with the 
free boundaries. In section V we present results for a vortex moving in diagonal JJA. Section 
VI contains a summary of our results and some conclusions. 



II. MODEL EQUATIONS 

In JJA vortices are represented by patterns of eddy current around a plaquette. We 
consider dc-biased JJA in the classical regime defined by Ej ^> E c = e 2 /2C, where Ej is 
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the Josephson coupling energy and E c the characteristic charging energy of a junction, e the 
electron charge, and C the mutual capacitance of the junctions. In this case, and in zero 
applied magnetic field, the individual junctions in the JJA can be modeled by the resistively 
and capacitively shunted junction (RCSJ) model [pOjl , defined by the total bond current 
i(r, r') between nearest neighbor sites r and r' as 

z(r, r') = pj(r, r') + 6(r, r') + sin[0(r, r')]. (1) 

Here the dots represent time derivatives. The three contributions to i(r, r') are the displace- 
ment, the dissipative and the superconducting currents, respectively. The phase difference 
across a junction is 9(r, r') = 9(r) — 9(r'). The currents are expressed in units of the junction 
critical current I c ; time is measured in units of the characteristic time l/u c = h/(2eR n I c ), 
(3 C = (uc/ujp) 2 is the Stewart-McCumber damping parameter P0| , with the plasma frequency 
u p defined as uo 2 p = 2eI c /%C, and R n is the junction's normal-state resistance. The RCSJ 
model has a limitation. It does not take into account that the resistance of a real Josephson 
junction with high (3 C is voltage-dependent: for voltages below the gap voltage the effective 
resistance is determined by the quasi-particle resistance, which can be orders of magnitude 
larger than the normal-state resistance [|1^ . 

In Fig. [I] we show the geometry of the square array to which most of the results discussed 
in this paper apply. Each superconducting island is connected to four neighbor islands via 
identical Josephson junctions. The bias current is fed in at the lower boundary and taken 
out at the upper boundary. 

When we neglect the self-induced magnetic field produced by the current flowing in the 
array (infinite magnetic penetration depth A), the array dynamics is described by Eq. (|TJ), 
together with Kirchhoff 's current conservation condition 

^z(r,r + a) = z ext (r), (2) 

a 

imposed at every island r (the summation is over all nearest neighbor islands r + a). We 
integrate the coupled equations (|l|) and (0) using a fast algorithm discussed in plf. 



When we include the self-induced fields, thus including screening effects, we redefine 
9(r,r') in ([I]) to be the gauge invariant phase difference across a junction: 9(r, r') = 9(r) — 
9(r') — 27T74(r,r'). The bond frustration variable A(r, r') is defined by the line integral of 
the vector potential A: 

where $o = h/2e is the elementary flux quantum. The vector potential is time-dependent, 
as is the total magnetic flux $(R, t) at plaquette R: 

<£A(r,x> t t). (4) 

Here V(R) denotes an anti-clockwise sum around the plaquette R. The time-dependent 
magnetic flux is induced by all the currents flowing in the array via the Faraday- Ampere 
laws. In zero external magnetic field, we may write 
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$(R,t) = £r(R,r,a)i(r,r + a,t). (5) 



r.a 



T is a matrix that explicitly depends on the geometry of the array and the junctions. The 
standard inductance matrix is then given by the discrete curl of F. As discussed in more 
detail in Ref . |22J , in the linearized approximation the strength of the self-induced fields is 
governed by the perpendicular magnetic penetration depth 



or equivalently by the ^-parameter 



a uj c ' 



(7) 



where uj<$> = R n /(fiQa), and fiQ is the vacuum magnetic permeability. We set the lattice 
constant a equal to unity and use the parameter A to indicate the degree of screening. 

In our calculations we use the full-range inductance matrix ("Model C" in Ref. |[22|| ). 
The inductance matrix elements are computed using the model approximations discussed in 
Ref. After choosing the temporal gauge, equations (|l|) to (|[) can be combined into one 
(vector) equation that governs the complete dynamics (see Ref. fl22fl). 

The vorticity n(R) of a plaquette R is defined as 

27m(R) = 2tt$(R) + £ (0(r) - 0(r') - 2nA(r, r')). (8) 
Here the gauge invariant phase difference 9(r) — 9(r') — 2irA(r, r') is taken between — 7r and 

+7T. 

We would like to model the vortex motion in the array by a coarse-grained equation 
of motion as in the Bardeen-Stephen model for flux flow in continuous superconductors. 
The dissipative currents in the junctions, modeled in Eq. ([I]) by currents through ohmic 
shunt resistors, give rise to a viscous force on a moving vortex. Furthermore, in capacitive 
arrays (/3 C > 0), the electromagnetic energy stored in the junction capacitors due to the 
vortex motion, can be interpreted as the kinetic energy of a massive vortex P 4T4]j23|1 . For an 



infinite array, and for infinite magnetic penetration depth, a coarse-grained model equation 
in terms of a single continuous vortex coordinate x reads: 

Mx + r/x = it, + id sin(27rx), (9) 

where M = ir(3 c and r\ = 7r for a square array P- |T4] , p3|1 . This equation describes the vortex 
as a point particle with mass M that, driven by a (Lorentz) force proportional to %, moves 
through a sinusoidal pinning potential and experiences a viscous damping force with constant 



viscosity coefficient r\. An estimate for i& in a square lattice gives id ~ 0.10 |24|. In a previous 
paper ||, we have deduced an alternative phenomenological vortex equation of motion with 
a velocity-dependent vortex viscosity. In contrast to the model with constant viscosity given 
in Eq. (Q), the nonlinear model gives a qualitatively correct account of the numerical results 
for the current- voltage characteristics calculated using the full set of equations (p. The 
equation reads 
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M(j3 c )x + - + ^ x = i b + i d sm(2nx). (10) 
Here the phenomenological constants A, B and M are found to be weakly dependent on j3 c . 



In the model equations @ or (10 ) the vortex dynamics is described in terms of a single 
particle with coordinate x. In writing these equations we assume that, as the vortex moves, 
it does not couple to other excitations in the array. As a result, it is assumed that the 
electromagnetic energy associated with the moving vortex, and interpreted as its kinetic 
energy, can only be absorbed by the viscous medium and not transmitted to other dissipative 
modes in the array. The limited validity of this assumption is apparent for example from 
numerical simulations nMfl, where no ballistic vortex motion was found when switching 



off the bias current. Furthermore, it was found in experiments [[IIJ and simulations that 
the vortex viscosity does not decrease as l/R n when increasing R n . The non-zero vortex 
viscosity in the underdamped limit was understood as due to the coupling of the moving 
vortex to charge oscillations ('spin waves') [pj|7 yi4]Jl6l . 

For a finite array, an extra position-dependent force due to the vortex-boundary inter- 
action should enter at the right-hand side of these equations: 

M'x + rjx = %b + id sin(27rx) — F(x), (11) 

and 

M((3 c )x + - + ^? c j x = i b + i d sin(27rx) - F(x). (12) 

To compute this force, one has to take into account all the infinitely many image vortices 
that are produced by the array boundaries. When we take the origin at the left-hand side 
boundary of the array, this leads to the following form of the force: 

F M = -;s ln <T sm <T»' < 13 > 



In order to avoid the singularities at the boundaries in (|lq) , we use a cutoff: for distances 
to the boundaries smaller than half a lattice constant, we set F(x) = 0. Note that this 
cutoff prescription fixes the maximum attractive force due to the boundary on the vortex 
at -fmax = F(x = \). In Ref. (TB| it is argued that the vortex mass increases close to a 



free boundary. It was found that this increase is only noticeable within one lattice constant 
from the boundary. We do not include this quantitative correction in the models (|TTD and 
( [12]) , as the arbitrariness of the cutoff prescription already makes that the results are only 
qualitative in nature. 

III. RESULTS 

A. A = oo 

In Fig. we present a summary of the results for a 16 x 16 square array with A = oo for 
different values of f3 c . At the bottom we show the results of the simulations for (3 C = 0, in 
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which the vortex motion is purely viscous, and the vortex mass is zero. We have considered 
all current values between if, = 0.0 and if, = 1.0, with a grid of Az = 0.01. For this system 
size the vortex depinning current is between 0.11 and 0.12, so for currents if, < 0.11 the 
vortex is pinned to the middle plaquette of the array, and for i b > 0.12 it is depinned (If 
one increases the system size the depinning current will reach the value id = 0.10 estimated 



in Ref. [24| for an infinite array). When the vortex is depinned and moves towards one of 



the free boundaries, the approach to its image on the other side of the boundary causes it 
to accelerate. When the vortex reaches the boundary, it leaves the array, or, equivalently, it 
is annihilated by its image antivortex. We denote this behavior as type A. 

Above i b = 0.96 the f3 c = dynamics looses its single-vortex character after a short 
transient time, when additional vortices are generated in other rows. In fact the measured 
voltage is due to the motion of vortices in all rows. As we want to focus on the single-vortex 
dynamics in this paper, we will not go further into this type of motion here. 

For nonzero f3 ci we started the simulations at t — increasing the current linearly from 
zero to if, in one RC time of the junctions. The relaxational oscillations due to this increase 
in the current decay while the vortex is moving towards the boundary. We checked that 
the behavior at the boundary is not seriously affected by these relaxational oscillations by 
comparing the results to those in a 32 x 16 array. 

For (3 C = 2, the results for the vortex propagation are similar to the (3 C = results. The 
only difference is that the regime with type A behavior now ends at i = 0.86 due to the onset 
of row switching, which means that the vortex switches one or more rows of longitudinal 
junctions (i.e. junctions in the current direction) into the resistive state jIH3l|2"5"|. 

For f3 c > 2.5 a current range opens up in which the vortex is reflected as an antivortex 



at the boundary in a way a soliton reflects in a long Josephson-j unction (LJJ) |T7,18 



Alternatively, one may say that the vortex and its antivortex image pass through each 



other, analogously to the non-destructive soliton-antisoliton collisions in a LJJ [|17i|l9[ . The 



antivortex in turn is reflected at the opposite free boundary as a vortex. This sequence 
repeats so that the vortex/antivortex never escapes from the array, thus producing a non- 
zero time-averaged voltage perpendicular to the vortex motion. We denote this as type 
B behavior. In a LJJ a similar type of soliton motion gives rise to the first zero-field 
step | 20fl . Although there are important differences between the properties of solitons in a 
continuous junction and the properties of vortices in the discrete array, it appears, in the 
regime considered here, that their reflection properties at a boundary as well as their collision 
properties are similar in many respects. For a two-dimensional array of Josephson junctions 
the non-destructive collision of a vortex-antivortex pair has been alluded to by Nakajima 
and Sawada, in the context of a model including the self-inductances of the plaquettes [f2(| . 

In Fig. ^| we show the voltage versus current characteristic for (3 C = 10. We note 
that for the calculation of this current-voltage characteristic we use the same initial phase 
configuration for all the bias currents considered. Above the vortex depinning current, we 
start the calculation of the time-averaged voltage only after a sufficiently long time interval 
for the vortex to reach the boundary. As a result, in the case of type A behavior we measure 
no voltage, as there is no dissipation after the vortex escapes from the array. In contrast, the 
type B regime has a non-zero voltage. The current- voltage characteristic has a considerable 
slope in this regime, in contrast to the first zero-field step in a LJJ, which is much flatter, 
due to the fact that the voltage is limited by the maximum soliton velocity. For even higher 
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currents, we enter the row-switching regime with a large voltage increase. 

We interpret the type B behavior as being the result of the inertia, or kinetic energy, 
carried by the vortex. The attractive interaction between the vortex and its image provides a 
potential well from which the reflected antivortex has to escape after the reflection/collision, 
in order to travel towards the opposite boundary. The Lorentz force on the antivortex due 
to the applied bias current is in itself not sufficient to pull it out of the well. In addition, 
the vortex needs to have a minimum kinetic energy in order to escape: this will be the 
case if both the vortex mass M = irj3 c as defined in equation (|9|) and the vortex velocity 
(monotonically increasing with are sufficiently large. In Fig. [| we show snapshots of the 
vortex configurations at different times for (3 C = 10 and if, = 0.49. In Fig. |5| we show the 
time-dependent voltage across the array indicating the times at which the snapshots shown 
in Fig. |] are taken. The large voltage fluctuations for short times are due to the additive 
contributions of the relaxational oscillations of all the individual longitudinal junctions in 
the array as the bias current is changed from zero to 0.49 in the time interval < t < 10. 
For t > 100 these oscillations have relaxed sufficiently and the vortex contribution to the 
voltage becomes dominant. As the vortex approaches the boundary its velocity increases 
and so does the voltage. At the reflection, the vortex velocity jumps from a maximal value 
(just before the reflection) to a minimal value (just after the reflection). 

To verify the interpretation of the vortex reflection as a non-destructive collision with 
the image antivortex, we have investigated the collision properties of a vortex-antivortex 
pair explicitly. This entails simulations of a 16 x 32 square array with periodic boundary 
conditions in the direction of the vortex motion, with one vortex and one antivortex separated 
by 16 lattice constants present in the initial phase configuration. This situation is the 
explicit realization of the image-vortex system corresponding to the 16 x 16 array with free 
boundaries and a single vortex. Below the row-switching threshold, we find that indeed the 
corresponding vortex dynamics, i.e. destructive collision A' and non-destructive collision 
B', takes place. The transition between the types A' and B' behavior is at a current value 
close to the one between types A and B in the finite array with one vortex, thus providing 
a phenomenological explanation of the results. 

For currents in the type A' regime close to B', the destructive collisions become a two- 
stage process: first the vortices collide non-destructively, drift three lattice constants apart 
and then fall back onto each other and annihilate. In the LJJ-soliton language this is called 
a decaying breather mode. In the finite array the corresponding process in the A range 
(near to the B range) consists of a reflected antivortex that falls back over the edge. We 
interpret this process as evidence for a non-zero vortex mass: the vortex moves up-hill (in 
the potential landscape) until it reaches a turning point where the kinetic energy is used up, 
and then falls back in the potential well. 

When j3 c increases, the lowest current resulting in type B behavior decreases, as shown 
in Fig. 0. The current for row-switching decreases as well, similar to the (3 C dependence 
of the row-switching threshold in previous simulations with periodic boundary conditions 
In the next section we will interpret the type A versus type B behavior in terms 
of the coarse-grained model equations given in Eqs. (|ll]) to (p~3|) of a vortex interacting 
logarithmically with its image(s). 

For (3 C > 50 we observe the sequence A^B->A->Basa function of bias current: 
the type B regime splits up into two pieces separated by a second type A regime. The 
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second type A regime is due to the presence of charge oscillations on the shunt capacitors, 
which become larger in amplitude for increasing f3 c . At the reflection, the vortex velocity 
becomes so large, that additional dissipative modes in the array are excited in the form of 
local charge oscillations on the shunt capacitors. These oscillations interfere with the motion 
of the antivortex, which as a result is slowed down and cannot escape from the potential 
well. For higher currents, there is a second B region. Here the antivortex has a higher 
velocity, and it is able to survive the coupling to the charge fluctuations. The magnitude 
of the second A regime grows with /3 C , and we have found that for (3 C = 500 there are no 
B regimes left. We note that the break up of the type B regime into two pieces and the 
disappearance of the reflection for very low damping have no counterpart in the context of 
a soliton moving in a LJJ. 

One might ask if the collision can still be non-destructive if the vortices collide along 
a direction that is not parallel to one of the coordinate axes. We have simulated a finite 
16 x 16 square array with f3 c = 10 containing a vortex and an antivortex that are inserted 
not in the same row but in adjacent rows. The results show a reduced but non- vanishing 
current regime (0.63 < % < 0.66) for reflections for both vortices at the edges plus a 'new' 
non-destructive collision in the middle of the array, during which the vortices interchange 
rows, as shown in Fig. |5|. 



B. Finite A regime 

Thus far we have studied the reflection properties of the vortex in an array with zero 
self-induced magnetic fields, or equivalently, with magnetic penetration depth A much larger 
than the array size. From the point of view of for example arrays made of niobium junctions, 
in which the self- induced magnetic fields are non-negligible |[27||, it is important to see how 



these fields influence the reflection properties of a vortex. We have therefore studied the 
vortex propagation in 16 x 16 square arrays for a range of A values at j3 c = 10 and 100. We 
find that the reflection still occurs for finite A. The results are shown in Figs. and |8]. As 



discussed by Phillips et al. ||27|| ) the vortex depinning current is enhanced by the self- induced 
fields. For example, for our 16 x 16 array with A = 0.3 and j3 c = 100 we find that the vortex 
depins for currents larger than or equal to ib = 0.38. For (3 C = 10 and in the type-II regime 
(A larger than a few lattice spacings) the width and the position of the B regime is somewhat 
insensitive to the A value as one can see in Fig. [F[ As we see in Fig. ||| for (3 C = 100, the 
second A regime shrinks with decreasing A, and disappears for A < 5. For A < 1 (type-I) 
we find that the B regime shrinks, both for (3 C = 10 and f3 c = 100. 



IV. COMPARISON WITH THE PHENOMENOLOGICAL MODEL EQUATIONS 

In this section we will make a qualitative comparison of the A = oo results presented in 
section III with the results for the vortex motion in a finite array based on Eqs. ( |TTD and 

(0)- 

Naively, in Eq. (|13D the parameter range for x is < x < L. However, when the vortex 
reflects at x = L, the positive image vortex of the reflected antivortex has a coordinate 
x > L, and we can interpret the oscillating state (type B behavior) as the motion of a 
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positive vortex through the force field F(x), periodically extended to the complete real 
axis. A behavior then corresponds to trapping of the vortex around x = L (modulo L). 
Starting from x = 8 at each current value considered, we calculated the current-voltage 
characteristic numerically from Eqs. (|TT|) and ([H|) for L = 15, using the cutoff prescription 
mentioned above. The results for (3 C = 10 are shown as the dashed line in Fig. [| The model 
underestimates the threshold for type B dynamics. As this threshold is very sensitive to the 
type of cutoff used, we can not draw quantitative conclusions about the vortex mass from 
this result. 

We have determined the threshold currents ia^b as a function of f3 c from the model 
equation flUD , as well as from the same equation with the nonlinear friction term from Eq. 



(|i2|). We find that the dependence of the A — > B threshold current on (3 C is qualitatively the 
same as found in the simulations. In the models the decrease of the threshold value is due 
to the increase in the mass parameter M = 2n/3 c , which increases the vortex kinetic energy. 
The larger this kinetic energy is, the easier it is for the vortex to escape from the effective 
potential well. 

We note that the way in which vortex inertia manifests itself, and hence the (possible) 
attribution of a mass, depends on the dynamical situation considered. In previous simu- 
lations []IIfl|3J], in contrast to this work, the vortex mass was probed by switching off the 
bias current. Changing the current gives rise to enhanced local relaxational oscillations in 
the junction phases near to the vortex center, or in other words, to a coupling to other 
dissipative modes in the array. Instead of allowing the vortex to continue its motion, the 
electromagnetic energy stored in the capacitors then leads to local oscillations. For large 
f3 c these force the vortex center to oscillate back and forth a couple of times between two 
adjacent plaquettes, as shown for (3 C = 2500 in Fig. 10 of Ref. 0. During the decay of these 
oscillations the electromagnetic energy is dissipated. The distance traveled by the vortex 
after the current has been turned off, is zero. Therefore the mass attributed to the vortex 
in this situation is zero or very small. Similarly, probing the vortex mass by looking at the 
hysteresis in the I-V characteristics [0,0, which also involves changes in the current bias, 
leads to zero or very small phenomenological mass M(f3 c ) in ([10|) for moderate values of f3 c 

On the basis of the model interpretation, one would expect that, for some currents just 
below the transition to type B behavior, the reflected antivortex is not falling back over 
the boundary, but it is re-trapped by the lattice. This re-trapping is then the effect of the 
combination of the pinning modulation of the effective potential and the energy loss due to 
friction. We have looked for such a re-trapping process in the simulations for (3 C = 100, and 
indeed found it for if, = 0.2462 (16 x 16 array). The reflected antivortex is re-trapped in the 
second plaquette from the boundary. 



V. VORTEX REFLECTION IN DIAGONAL ARRAYS 

An important question is how robust the vortex reflection at a boundary (or the non- 
destructive collision of a vortex-antivortex pair) is for other array geometries. We have 
studied the case in which the current is injected not along one of the coordinate axes, but 
along the — x and y direction. See Fig. |9|, that shows a 15 x 15 diagonal array, in which 
the initial position of the vortex is indicated by an open circle. The vortex will move at 
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45 degrees to the boundary, due to the diagonal current bias, which allows us to study the 
reflection behavior for this angle of incidence. We also find soliton-like reflection in this case, 
albeit the width of the current regime is smaller, and vanishes again already at a (3 C value 
somewhere between 50 and 100. Although the vortex and its image collide at an angle of 
90 degrees, the reflecting antivortex still retraces the path of the incident vortex. This is 
basically because the antivortex has to travel at right angles to the applied bias current. 
An interesting question would then be what happens for a reflection at zero current bias. 
To study this within the setting of our model simulation, one would need a nonzero current 
first to depin the vortex and then switch the current off just before the reflection. However, 
it is known |],0|4[] that switching off the bias current is not followed by a propagation of the 
vortex for /3 C in the range in which we find vortex reflection. Therefore, switching off the 
bias current just before the reflection does not seem to be an option to study reflection at 
zero bias current. 



VI. SUMMARY AND DISCUSSION 

In conclusion, we have performed RCS J-model simulations of a finite Josephson-junction 
array containing a single vortex in a dc current bias. For moderate damping, we find current 
regimes where the vortex reflects at the boundary, moving back as an antivortex towards 
the opposite boundary, where the dual reflection process takes place, leading to a stationary 
oscillatory state and a non-zero time-averaged voltage across the array. This oscillatory 
motion can be viewed as a discrete analogue of the bouncing of a soliton in an underdamped 
long Josephson junction, which gives rise to the first zero-field step ||17|| . The long Josephson 



junction or Josephson transmission line can be used as a vortex flow transistor or oscillator in 
many different devices |[28|| . Recently, Van der Zant and Orlando p9| explored the possibility 
of using a discrete one-dimensional parallel array of underdamped junctions as a vortex flow 
transistor. Some time ago it has been shown numerically [|30| as well as experimentally 



31,32] that vortices in a ID parallel Josephson array have soliton-like collision properties. 
Our results show that the soliton-like collision and reflection also occur in two-dimensional 
arrays, and that the first zero-field step has a discrete analogue in JJA, which suggests a 
possible use of the two-dimensional array as a vortex-flow device. An important difference 
found with respect to the soliton motion in LJJ for very low damping, is the fact, that in 



the array the reflection is absent for values of (3 C of the order of 500 and larger p9 |. 

When considering a vortex-antivortex pair in an array with periodic conditions, the vor- 
tices collide non-destructively for appropriate bias currents. This is the behavior equivalent 
to the reflection of one vortex at a free boundary. Our numerical results can be interpreted 
in terms of a macroscopic model equation for a massive vortex. We also studied the vortex 
propagation in a diagonal array where the vortex has an angle of incidence of 45 degrees to 
the boundary, and found that, for appropriate bias currents, reflection takes place here as 
well. Finally, we also studied the effect of self-induced magnetic fields, using a model that 
takes into account the full inductance matrix of the array and found that the same type of 
vortex reflection regimes as in the A = oo case are present for a range of values of A and f3 c . 



10 



ACKNOWLEDGMENTS 



We thank J.E. Mooij for suggesting the topic of this work. We also thank him and W. 
Elion, P. Hadley, A. van Oudenaarden and H. van der Zant for discussions. This work was 
supported in part by the Dutch organization for fundamental research (FOM). The work of 
JVJ has been partially supported by NSF grant No. DMR-9521845. 

* Present address: Institut fur Theoretische Physik, Universitat Wiirzburg, Am Hubland, 
97074 Wiirzburg, Germany. 



11 



REFERENCES 



[1] P. A. Bobbert, Phys. Rev. B 45, 7540 (1992). 

[2] U. Geigenmuller, C.J. Lobb, and C.B. Whan, Phys. Rev. B 47, 348 (1993). 
[3] W. Yu, K.H. Lee, and D. Stroud, Phys. Rev. B 47, 5906 (1993). 
[4] W. Yu and D. Stroud, Phys. Rev. B 49, 6174 (1994). 

[5] T.J. Hagenaars, P.H.E. Tiesinga, J.E. van Himbergen, and J.V. Jose, Phys. Rev. B 50, 
1143 (1994). 

[6] T.J. Hagenaars, P.H.E. Tiesinga, J.E. van Himbergen, and J.V. Jose, Quantum Dy- 
namics of Sub micron Structures, Vol 291 of NATO Advanced Studies Institute, Series 
E: Applied Sciences, edited by H.A. Cerdeira et al (Kluwer, Dordrecht 1995), p. 617. 

[7] U. Geigenmuller, preprint (1995). 

[8] E. Simanek, Inhomogeneous Superconductors, Oxford University Press, New York 1994. 
[9] A.I. Larkin, Yu. N. Ovchinnikov, and A. Schmid, Physica 152B, 266 (1988). 
[10] U. Eckern and A. Schmid, Phys. Rev. B 39, 6441 (1989). 

[11] U. Eckern, in Applications of Statistical and Field Theory Methods to Condensed Matter, 
Vol. 218 of NATO Advanced Studies Institute, Series B: Physics, edited by D. Baeriswyl 
et al. (Plenum, New York 1990), p. 311. 

[12] M.S. Rzchowski, S.P Benz, M. Tinkham, and C.J. Lobb, Phys. Rev. B 42, 2041 (1990). 

[13] T.P Orlando, J.E. Mooij, and H.S.J, van der Zant, Phys. Rev. B 43, 10218 (1991). 

[14] U. Eckern and E.B. Sonin, Phys. Rev. B 47, 505 (1993). 

[15] H.S.J, van der Zant, F.C. Fritschy, T.P. Orlando, and J.E. Mooij, Europhys. Lett. 18, 
343 (1992). 

[16] H.S.J, van der Zant, F.C. Fritschy, T.P. Orlando, and J.E. Mooij, Phys. Rev. Lett. 66, 

2531 (1991); Phys. Rev. B 47, 295 (1993). 
[17] T.A. Fulton and R.C. Dynes, Solid State Comm. 12, 57 (1973). 
[18] K. Nakajima, Y. Onodera, and Y. Ogawa, J. Appl. Phys. 47, 1620 (1976). 
[19] D.W. McLaughlin and A.C. Scott, Appl. Phys. Lett. 30, 545 (1977). 
[20] see for example: K.K. Likharev, Dynamics of Josephson Junctions and Circuits (Gordon 

and Breach, New York, 1986). 
[21] H. Eikmans and J.E. van Himbergen, Phys. Rev. B 41, 8927 (1990); A similar but even 

faster algorithm was introduced in D. Dominguez, J.V. Jose, A. Karma, and C. Wiecko, 

Phys. Rev. Lett. 67, 2367 (1991). 
[22] D. Dominguez and J.V. Jose, Int. J. Mod. Phys. B 8, 3749 (1994). 
[23] H. Eikmans and J.E. van Himbergen, Phys. Rev. B 45, 10597 (1992). 
[24] C.J. Lobb, D.W. Abraham, and M. Tinkham, Phys. Rev. B 27, 150 (1983). 
[25] H.S.J, van der Zant, F.C. Fritschy, T.P. Orlando, and J.E. Mooij, Physica (Amsterdam) 

165&166B, 969 (1990); and H.S.J, van der Zant, C.J. Muller, L.J. Geerligs, C.J.P.M. 

Harmans, and J.E. Mooij, Phys. Rev. B. 38, 5154 (1988). 
[26] K. Nakajima and Y. Sawada, J. Appl. Phys. 52, 5732 (1981). 

[27] J. R. Phillips, H. S. J. van der Zant, J. White, and T. P. Orlando, Phys. Rev. B 47, 
5219 (1993). 

[28] N.F. Pedersen, IEEE Trans. Magn. 27, 3328 (1991). 

[29] H.S.J, van der Zant and T.P. Orlando, J. Appl. Phys. 76, 7606 (1994). 

[30] K. Nakajima and Y. Onodera, J. Appl. Phys. 49, 2958 (1978). 

[31] A. Fujimaki, K. Nakajima, and Y. Sawada, Phys. Rev. Lett. 59, 2895 (1987). 



12 



[32] K. Nakajima, H. Mizusawa, Y. Sawada, H. Akoh and S. Takada, Phys. Rev. Lett. 65, 
1667 (1990). 



13 



FIGURES 



FIG. 1. Square array geometry used in the simulations, illustrated with a 8 x 8 array. Junctions 
are denoted as crossed bonds. Free boundary conditions are imposed in both directions, while the 
current bias is applied along the y-direction. 

FIG. 2. Results for different values of (3 C for a 16x16 square array with free boundaries. The 
dashed lines represent the current ranges for which type A behavior (vortex escapes from array) 
was found. The thick lines denote the type B ranges (trapped vortex oscillating in the array) and 
the thin full lines the row-switching regimes. 

FIG. 3. Current-voltage characteristics for /3 C = 10, from simulations of a 16 x 16 square array 
with one vortex (full line). (V)t is the time-averaged voltage across the array in the current 
direction, in units of I c Rn an d normalized by the number of longitudinal junctions. The dashed 
line is the I-V characteristic computed from the linear viscosity model (11). 

FIG. 4. Snapshots of the vorticity distribution Eq. (8) in a 16 x 16 square array with (3 C = 10 
and % = 0.49, showing type B dynamics. In the frame labeled with the vortex (black square = 
vortex center plaquette) is moving towards the right free boundary. After reflection it travels as 
an antivortex (white square) towards the opposite boundary (frames 1 and 2), where it is reflected 
again as a positive vortex (frame 3). The snapshots correspond to the instants indicated in Fig. 5. 

FIG. 5. Voltage, normalized as in Fig. 3, versus time across a 16 x 16 square array with one 
vortex, for (3 C = 10 and i = 0.49. The times labeled with to 3 correspond to the respective 
snapshots shown in Fig. 4. 

FIG. 6. Snapshots of the vorticity distribution for a 16 x 16 array with free boundaries and 
% = 0.65, with one positive vortex and one antivortex moving in adjacent rows. The notation is as 
in Fig. 4. Both vortices reflect at the edges (between frames and 1), and collide constructively 
in the middle of the array (between 2 and 3), interchanging rows. 

FIG. 7. Results for different values of A for a 16x16 square array with free boundaries at 
(3 C = 10. The notation is the same as in Fig. 2. 

FIG. 8. Results for different values of A for a 16x16 square array with free boundaries at 
(3 C = 100. The notation is the same as in Fig. 2. 

FIG. 9. Diagonal array geometry used in the simulations, illustrated with a 15 x 15 array. 
For clarity we omitted the crosses on the bonds. In both directions free boundary conditions are 
imposed, while the current bias is applied along a diagonal direction. The open circle denotes the 
initial position of the vortex. 
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